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The literature in astrodynamics abounds with •perturbation techniques 
for satellite orbits. Various formulations have been generated in terms of 
orbit elements, the satellite position and velocity vectors, or combinatiotis 
thereof. The computational effectiveness of any perturbation scheme depends 
largely on the definitions used for the dynamic state variables. Some meth- 
ods are aimed at long-range predictions and orbit lifetime studies, others 
at short-range predictions for guidance. This paper may serve as an intro- 
duction to this field for the nonspecialist, in that it reviews the classical 
variation-of-parameters technique and discusses several engineering analy- 
ses that were generated in the post-Sputnik era. It also points to some con- 
nections between these relatively simple approaches and more elaborate 
methods of celestial mechanics. Thus it may contribute toward a comparison 
of several ''professional" approaches whose relative merits are often de- 
bated among experts. 

I. INTRODUCTION 

This paper is a discussion of various perturbation techniques for 
satellite orbits which were investigated by the author and his colleagues 
during the past few years. The effort began with a tutorial "orbit 
seminar" several years ago and it seemed appropriate to collect some 
of this material here as a companion paper for R. B. Blackman's "Meth- 
ods of Orbit Refinement." 

It is a symptom of our times that aerospace engineers are taking a 
new look at the established methods of dynamical astronomy. The 
orbital geometries and vehicle characteristics encountered with arti- 
ficial celestial bodies often require departures from the formulations of 
classical astronomy and, in fact, have stimulated several new (or at 
least independent) approaches during the post-Sputnik era. The number 
of publications in this time has been formidable, and in many discus- 
sions the names attached to various formulations serve as passwords 
for the ideas they represent. The uninitiated find themselves at a loss 
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concerning the methods that stand behind these names, their degree of 
originality, and their relations with each other. 

In view of this situation the following article is addressed to two kinds 
of readers: 

(i) The newcomers in the field of orbital mechanics who seek a tutorial 
survey and an introduction to some of the literature. A bare minimum 
of definitions is given for their benefit; a discussion of basic order-of- 
magnitude relations and certain intuitive notions which would strengthen 
the beginner's grasp of the physical problem had to be omitted for lack 
of space but can be found in the literature. 1,2 

(it) The specialists in orbital mechanics who have not had occasion 
to correlate some of the better-known contributions in the literature 
and who may find this work a step in that direction. Typical issues 
in such comparisons are the choice of coordinates, the accuracy and 
elegance achieved by various transformations of the variables, and the 
precision obtainable from series expansions of the solution in terms of 
various small parameters. 

The simultaneous need for conciseness of presentation and discussion 
of certain analytic detail presents somewhat of a dilemma. As a com- 
promise, much of the development between the explicitly quoted re- 
sults is covered in a descriptive way and the reader is referred to the 
literature for all standard derivations. 1 " 4 Most of our discussion con- 
cerns orbits of moderate eccentricity, which are representative of satel- 
lite missions. However, in many places an extension to the highly ec- 
centric orbits of space probes follows readily. 

We begin by devoting Section II to a statement of the fundamental 
equations of motion, the definition of so-called orbit parameters, and a 
description of various disturbing functions. Section III summarizes 
the classical treatment of satellite perturbations as gradual changes 
of the orbit parameters. (From a general point of view, this formulation, 
due to Lagrange, is derivable from the canonical systems governing the 
satellite problem.) It is hoped that this covers a sufficient amount of 
standard material to introduce the concepts and the parlance of orbital 
mechanics. 

In Section IV we examine several perturbation methods for aero- 
space applications which are based on variously defined spherical and 
moving Cartesian coordinates. This includes the well-known contri- 
butions by Blitzer et al., Anthony et al., and Roberson. They could 
serve as an introduction to the discussion of more elaborate formula- 
tions by King-Hele et al. and Brenner et al. In Section V we treat one 
more formulation in this general category which was specifically de- 
signed for guidance studies. 
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A logical continuation of this paper would cover the methods of 
Breakwell et al. and Diliberto, Kyner's averaging technique, and the 
one suggested by Struble. Ultimately the hierarchy of perturbation 
methods leads to the Hamilton-Jacobi techniques expounded by Brou- 
wer, Garfinkel, and Vinti. These represent a very popular approach to 
higher-order perturbations and the coupling between simultaneous 
disturbances of satellite orbits. 

II. PRELIMINARIES AND DEFINITIONS 

We remember that the underlying phenomenon of undisturbed satel- 
lite motion (in a central force field, i.e. around a spherically symmetric 
body) is Newton's law of inverse square attraction. In a Cartesian co- 
ordinate system this spells out to be 

x = - (Gmrx/r 3 ) (1) 

y = - (Gm,.y/r 3 ) (2) 

z = - (Gm c z/r 3 ) (3) 

where G is the universal gravitational constant, m e the central mass, 
and we shall usually take Gm v = k for brevity. m e shall be the mass of 
the earth in all our discussions. [Strictly speaking, formulas (1) to (3) 
should show the sum of m c and the satellite mass instead of just m r .] 
/• is the distance from the origin, and dots indicate time derivatives. 
The .{•-// plane is usually taken to coincide with the equator, while the 
positive x axis points to the vernal equinox. The above equations 
simply state that each acceleration component is due to the correspond- 
ing component of the gravitational attraction — the minus signs indi- 
cating a direction toward the origin. The solutions of (1) to (3) are 
the well-known Kepler orbits — ellipses, parabolas, and hyperbolas. 

Such orbits can be conveniently described by a set of six parameters 
that give the plane of the motion, the shape of the orbit and its orien- 
tation in that plane, and the timing of the satellite motion along this 
path. These quantities may be considered the constants of integration 
for a solution of (1) to (3). A standard set of such orbit elements for 
elliptic motion is illustrated in Fig. 1. They are: 

a, the semimajor axis, 

e, the eccentricity, 

to, the argument of perigee, (4) 

i, the inclination, 

12, the nodal angle, and 

t, the time of perigee passage. 
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Fig. 1 — Standard set of orbit elements for elliptic motion. 



The last quantity establishes a time scale for the entire motion in that 
it serves as "epoch" and fixes one particular passage through perigee. 
If the satellite has swept out the angle / since that passage, the elapsed 
time is given by 



t - T = (a/k) 1 (2 tan 



u\i 



-[fe)M] 

- e (1 - C 2 ) J 



sin/ 



1 + e cos /jo 



(5) 



where t is the time pertaining to the position 0'. / is known as the true 
anomaly and (5) holds for all values of this angle. If we set/ = 2tt 
this corresponds of course to a full revolution around the orbit, and the 
elapsed time interval is 



T = 2tt (a/ky, 



(6) 



which is known as the anomalistic period. The instantaneous position 
0' can also be denned in terms of other angles, the so-called eccentric 
anomaly E or the mean anomaly M, which will be defined later. They 
can be related to time in similar ways. 

The parameters in (4) represent a typical set of orbit elements. The 
position and velocity vectors at some epoch are an alternative suggested 
by (l)-(3). Most astrodynamical theories use variations and combina- 
tions of all these, but from the general standpoint of analytical dynamics 
most sets of six parameters (if they are independent of each other) may 
be regarded as sets of canonic variables. Before we proceed to detailed 
formulations we examine briefly the various physical disturbances which 
cause these parameters to change in time. 
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2.1 The F'Jffect of Extraterrestrial Gravitation 

If we consider the attractions from masses other than the earth we 
speak of "extraterrestrial" gravitation. In the presence of a disturbing 
body I*, ( 1 ) becomes 

Giii r x n (x — Xp X p \ . . 

where 

m p = the mass of R, 

r ip - [(.c - Xpf + (y - ij„) 2 + (z - z p f]\ 
the distance from the satellite to R, and 

r„ = [Xp + Up + z P 2 ]\ 

the geocentric distance of I'. 

The corresponding equations for y and z are obvious. Now it is often 
convenient in analytical dynamics to express the disturbing terms in 
x, y, z as partial derivatives (d/7/d.r), (dR/dy), (dR/dz) of a disturbing 
function R. For the present case we would have 



s _ /',., f l xx p + yyp + ^1 



(8) 



as can be easily verified. 

We see that the ratio of the second term in (7) to the first is of the 
order m p r /m r r p ' — k. Typical values for k in the planetary system are 
equal to or less than lfP\ Its snuillness is vital to the entire rationale of 
a perturbation technique. 

2.2 The Effect of the Earth's Oblateness 

The potential field for the; nonspherical earth can be represented to 
various levels of accuracy by a series of spherical harmonics. If we 
restrict ourselves to terms with rotational symmetry about the polar 
axis, we obtain the following disturbing function 

& Gm e R \J , o • 2 \ 
R = - (1 - 3 sin <p) 

(9) 



+ — — (o sin ip — 5 sin <p) -\- 



where 
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R = the earth's equatorial radius, 
<p = the geocentric latitude of the satellite, 
,/ = 1.6239 X 10" 3 , and 
H = 6.04 X 10~ 6 . 

This two-term series is sufficiently accurate for our purposes. 

2.3 The Effect of Atmospheric Drag 

The resistance encountered by a satellite from the atmosphere is a 
subject of considerable uncertainty and continued research. For one 
thing, the density of atmospheric gases as a function of geographic 
location, altitude and time is not well known; moreover, the laws of 
interaction between a satellite and this rarefied medium are incompletely 
understood. Doubts exist as to the transition from a continuum be- 
havior of the atmosphere to the gas-kinetic regime and the extent to 
which electric interactions play a role. Nevertheless, the classical drag 
law yields useful results in many cases and we shall concentrate on it. 
We let 



F D = -(C D A/2)pv a 2 (10) 



where 



F D = the total drag force on the satellite 

A = the frontal area of the satellite 
C D = the drag coefficient 

p = the atmospheric density 

v a = the satellite velocity relative to the atmosphere. 

The monotonic decay of p with altitude covers approximately ten 
orders of magnitude within typical satellite altitudes and remains the 
subject of extensive study. The relative velocity v a is simply the difference 
between the satellite's inertial velocity vector v(x,y,z) and the rotational 
velocity of the atmosphere V = ra cos <p, where a can usually be taken 
as the earth's angular motion (the diurnal rate) and V always points due 
east. One is frequently justified in employing an approximate vector rep- 
resentation of (10) : 

F„^ (CbA/2) p »(V - v). (11) 

For typical earth satellites this force is at least two orders of magnitude 
smaller than the central attraction, i.e., k ^ 10 . 
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2.4 The Effect of Radiation Pressure 

As the reader knows, solar illumination exerts some pressure on 
every satellite. The magnitude of this force depends on the reflectivity 
and geometry of the satellite and, strictly speaking, on the distance 
from the satellite to the sun. It frequently suffices to represent this dis- 
turbance as a constant force /3 per unit mass and to note that it is many 
orders of magnitude smaller than the central gravity force. 

III. PERTURBATIONS IN THE ELEMENTS 

The six orbit elements [see (4)] were constants for the case of central 
inverse-square attraction. However, if any additional forces act on the 
satellite these parameters will be subject to change. To emphasize their 
time dependence we might write them as a(t), e(t) etc. In fact, their 
numerical values at any time t describe the ellipse the orbiting body 
would follow if all perturbations vanished as of that instant. This 
trajectory is obviously tangent to the actual flight path at t and is 
known as the "osculating" orbit. The relation between the satellite 
position in the osculating orbit and in the x, y, z frame follows from the 
geometry of conic section trajectories: 

a(l — e ) r . , . , ., 

x = — . /i cos / + h sin / 

1 4- c cos / 

U = ^~r — : } t'" 1 G0S / + nh sin /] ( 12 ) 

X *"i € COS 7 

a(l — e 2 ) r . , 

z = - — ■ Wi cos / + n 2 sin / 

1 + c cos / 

where l\ , k , ■ ■ • n-> are functions of i, w, and 12. Hence .r, y, z are repre- 
sentable as functions of a, e, i, oj, SI, /. [As mentioned with (5), we 
could also work in terms of the independent variable E or M instead 
of/.] However, the complete definition of the osculating orbit also entails 
that 

dx , . „ ,x df 

X = df ( <w ' co ' n '-" f t 

i, = | (<w>a/) | (13) 

dz t ■ ^ ^n df 
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where a, c, i, co, Q,, r are treated as constants. In other words, the ve- 
locity as well as the position in the osculating orbit are representative 
of the actual motion. This is the full extent of the "condition of oscula- 
tion." 

A large part of classical celestial mechanics has been based on the 
concept of osculating orbits, and during the post-Sputnik era Lagrange's 
classical treatment of the perturbations in these elements has been 
exploited ad ultimo. Its inclusion in this article is justified mainly by the 
need for completeness in an introductory survey such as this. It also 
serves as a point of reference for the nonclassical "perturbations in the 
coordinates" in the next section and for the Hamilton-Jacobi techniques 
frequently used by astronomers. 

In essence, the Lagrange method consists of transforming the basic 
equations 

GmeX . BR 

x = ~~ — r~ + t~ 
r dx 

Gm r %i . dR /- , x 

r 2 by 

Gm c z . dR 
r 3 dz 

to six first-order equations in the orbit elements and approximating 
their solutions by quadratures. Remembering that these parameters 
represented the constants of integration for the Kepler problem, ( 1 )- 
(3), we note that the transition to a{t), e(t), • • ■ is nothing but La- 
grange's "variation of constants" designed to accommodate the terms 
\dR/d(x } y,z)\ in (14). In the process of transforming (14) by means of 
(12) we avoid the occurrence of second derivatives of the orbit elements 
by demanding that 

dx dx . 
— = — , i.e., 
dt dt ' ' 

dx . . dx . . dx di . dx dx A , dx . n ( , -v 

da de di dt da dSl dr 

etc. This of course results in the reduction of the system of three second- 
order equations (14) to six first-order equations. Equations (15) are 
simply a restatement of (13), the condition of osculation. 

In principle, (14) could be transformed by a straightforward substi- 
tution of (12). In order to simplify the algebra, however, one may 
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symmetrize those equations to a form where all the labor reduces to the 
evaluation of such quantities as 

/ dx_ d± _ dx dx\ /dy_ d# _ dy_ dy\ 
\dctidaj detjdeti/ \doti datj daij dan/ 

( dz dz dz dz\ , , 

+ 1 — t— - t- r— ) = lai t, aj], 
\<>ai aaj aXj dan/ 

with i, j = 1 • • • 6. 

The shorthand symbol that we have adopted for this expression is 
known as a Lagrange bracket, and a, and aj stand for any two of the 
orbit elements. These brackets have the properties 

[«,-,«,] = 0, [a,-,*,-] = - [«>,«,] (17) 

and 

- [a, , a,] - 0, 

which make them useful devices in numerous manipulations of analytic 
dynamics. With their help the equations (14) become 

(18) 
(19) 

(20) 



2a dR 
a = ~ T dV 




o(l - e) dR 1/1- eVdR 
ke dr e \ ka / dco 


di 1 


r . dR _ dR 

die dil 


dt [fca(l - e 2 )]* sin* 



it- ,. ,., 1 ,,,. ■ M (2D 

[ka(l — e-)] 1 sm % dt 



i - e-y 1 



ka ) ( 



dR _ e cot i dR 
de 1 — e- di 



(221 



and a corresponding equation for f which will be discussed a little later. 
The five equations given here describe the changing geometry of the 
satellite orbit. All six equations together are known as Lagrange's 
planetary or "variational" equations. 

In (18)-(22) we assumed that the perturbing forces were conserva- 
tive, i.e., expressible as (dR/dx), (dR/dy), and (dR/dz). In some 
situations, as for example in the case of drag, this is not so. Under these 
conditions it is convenient to represent the disturbing force by com- 
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ponents S, T, N which are in the radial direction, the direction of in- 
creasing true anomaly, and normal to the orbit plane, respectively. The 
derivation leading to (18)-(22) can be repeated with the appropriate 
modifications to yield a set of differential equations with S, T, N in the 
right-hand sides. For example 



a = 2 



3 

a 



Ml - *)_ 



,S e sin/ + T (1 + e cos/)], etc. (23) 



These are known as Gauss's form of the planetary equations. We note 
that they contain the true anomaly as an independent variable. More 
will be said about this presently. 

It is possible to show that the planetary equations, especially in the 
last form, lend themselves to an alternative derivation which appeals 
to intuition. It is based on the idea that any continuously acting per- 
turbation may be interpreted as a sequence of infinitesimal impulses 
whose cumulative time response can be represented by a convolution 
integral. This approach leads to equations like (18) -(23) without the 
manipulations involving Lagrange brackets. 6 

Inspection of (18)-(22) shows that their nonlinear right-hand sides 
preclude an exact solution except for very special forms of R. Unfor- 
tunately, none of the perturbations encountered in nature fall into this 
category. One therefore resorts to a process of successive approxima- 
tions. 

Assuming that all disturbances represented by R are small in relation 
to the central attraction (i.e., 

d(x,y,z)J/ r m 

as discussed in Section II) we consider the solution for the undisturbed 
motion, i.e. the Kepler problem, as a "zero-order" approximation to the 
actual case. Let its parameters be denoted a (0) , e (0) •■•. If we insert 
them into the right-hand sides of (18)-(22), these equations reduce to 
quadratures yielding a first approximation to the effects of R on the 
orbit. These results are denoted a (1) , e (1) • • • and known as "first-order" 
perturbations. We observe that (a (1) /a (0) ), (e (1) /e (0) ), • • • = 0(k). In prin- 
ciple, this process can be repeated indefinitely by substituting a • • • 
into the right-hand sides to integrate for a (n) • ■ ■ . The limit is usually 
reached when the results for a, e • • • have settled or human endurance 
is exhausted. (The latter constraint may be eventually eliminated by 
computer routines for symbol manipulations.) The convergence of this 
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process to the exact solution has been established by PoincarS and is 
of fundamental interest to the mathematician. Suffice it here to say 
that the "smallness" of perturbations discussed in Section II should be 
such as to justify the iterative process. 

When the right-hand sides of (18)-(22) are written out explicitly 
for any particular case, they tend to become awkward because a trans- 
cendental angle-time relation like (5) enters. Since the geometric de- 
scription of a perturbation It (or S,T,N) usually involves an angle like 
one of the anomalies very directly, it is convenient to use one of them 
as independent variable. The time relation which interconnects the 
anomalies for an osculating orbit (Kepler's equation) can be stated in 
terms of the eccentric anomaly E, the true anomaly, / and the mean 
anomaly 71/ as follows: 

e(l — e 2 ) sin/ 



E — e sin E = 2 tan 



lo^y-a 



1 + e cos / (24) 

= (lc/a)'(t - r) = M. 

This may serve as a definition of E and M. The quantity {k/a)' = n 
is referred to as the "mean angular rate." 

If we work in terms of the true anomaly, we can write the left-hand 
sides of Lagrange's equations as 

da df 

" = if it • etc - 

where an expression must now be found for/. If we choose to consider/ 
as the true anomaly in some osculating orbit valid at time / , then it 
can be related to time by ( 24 ) in terms of the unperturbed elements 
a"" = a , e"" = e u etc. We call it an "unperturbed" anomaly and 
designate it by/' '. From (24) one finds 



*«» 



_Oo (1 — r (p ")' t 
which transforms ( 18)-(22) to 

J (1) r 7/, 2\3-li 

da a () (1 — c„ ) 



'Lf 






(l + e cos/ (O y (25) 



1 + c cos/ ') -2 ^,etc. (26) 



When integrated, these expressions represent first-order perturbations 
in terms of the unperturbed anomaly, i.e., a ( "(/ (0) ), etc. 

They suggest the following procedure for generating a first-order 
satellite ephemeris: if / is the starting epoch we evaluate a (l) (/ (0) ) etc. 
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between the limits / <0 ' and .A ". The time fa at the upper limit follows 
from (24) in terms of a , e • ■ • . Using fa and &i = a + a (/i ), 
etc., in (24), we find/! , the true anomaly for the new osculating orbit. 
Changing the notation from a { to ai , etc., and Ji to /i (0) , we can now 
repeat the procedure for the next integration interval. The "updating" 
of orbit elements in the right-hand side of (20) amounts to a partial 
allowance for higher-order perturbations, while the recalculation of / 
at the beginning of each step represents essentially a first-order pertur- 
bation of the true anomaly. 

Instead of doing the latter by discrete increments, we can work with 
a "perturbed" anomaly by differentiating (24) with proper allowance 
for the time dependence of a, e and t. Using (25) one finds that 

df_ f 2(1 - ef-a (SR\ ad ~ W gg 

rff«» j® z fc(l + e cos/) 2 \da) efc(l + e cos/) 2 de 



(27) 



+ n ^ 7 r l [d + cos/) tan'- 
fc*(l 4- e cos/) 3 L * 

/„ 2 1N • , (1 - e 2 )esin/cos/"| 
- (2e - 1) sin/ - 1+ecos/ j 

[ a(l - e 2 ) dR 1 (1 - A> dRl 

L fce 3t e\ oft /fluj 

where (dR/da) means that ^ is to be differentiated with respect to the 
semimajor axis wherever the latter appears explicitly but not when it 
is contained in n. This avoids the occurrence of a term with (t — t). 
An expression analogous to (27) can be derived in terms of S,T,N; see 
for example Ref. 7, p. 4. Now, it is immaterial in a first-order approxi- 
mation such as (27) whether we consider/ or / (0) as the independent 
variable in the right-hand side. Let us assume the former and use the 
symbol df/df w = p(f). Then (26) becomes 

d -C = — K (1 - e ° 2) 7 (1 + eo cos/)" 2 f , etc. (28) 

The integration procedure now runs between consecutive limits / , 
/i , • ■ • fj , w i tn updating being required only in the orbit elements. 
The corresponding epochs t, are of course computable by substituting 
Qj , «/ • • ■ fj into ( 24 )- Tne relative advantages of integrating the per- 
turbative equations in terms of / <0) or/ depend on the problem at hand. 
As we shall see later, the choice between an unperturbed or perturbed 
independent variable is available in most perturbation methods. 

Up to this point we have restricted our discussion to the first five 
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orbit parameters, which describe the geometry of the osculating orbit. 
Wherever t appeared, as in (24), we assumed that it would be avail- 
able from a suitable sixth equation. This parameter is needed to corre- 
late the independent variable, such as /, with time. 

An equation for r, corresponding to (18)-(22), can be obtained by 
the process outlined before, which yields 

2d 1 OR . ail - e 1 ) dR 



I: da ke de 



(29) 



Sometimes it is convenient to work with the slightly different parameter 
X = —tot. The differential equation for it reads 

„, x i dR 1 — e dR , v 

X = -2{a/ky — - — . (.30) 

da e(aA-) 5 de 

[A superficial comparison of these two equations gives the startling 
impression that (30) is obtained by multiplying (29) with —n, thus 
neglecting the n term that should appear. This term is really absorbed 
in the difference between (dR/da) T and (dR/da) x , as implied by the 
two equations; i.e. partial derivatives of R with respect to the scmimajor 
axis, holding t or x constant as required.] 

One could transform (29) and (30) to / (or E) as independent 
variable and pursue the quadrature as we did before. However, since 
d/7/da involves df/da (or dE/da), we notice from (24) that this in- 
troduces the factors tan -1 {[(1 - e)/(l + e)\ h tan (//2)j (or E) 
into the integrands for r (l) (or x'")- They can be quite awkward. 

Several devices have been developed to circumvent this difficulty. 
According to one approach we transform (29) from t to .1/ = n(l — r). 
The necessary compensating factors arise thereby which eliminate all 
aperiodic terms. Transforming the integrated equation back to r, we 
have 



On (1 — en ) /" 

1 / 

k J .' /o (o) 

(1 + eo cos/ ) ( - r— 5 a x + n T\ (31) 

2akf eo(A-ar») J dCo 



+ 9 ^«'(S)}*" 



where 

/7i = ( k/fti) , Wo = (k/a<j) , i\ — to + Ti , ai = a + cti , 

and (dft/dao) has the previously established meaning. Equation (31) can 



860 THE BELL SYSTEM TECHNICAL JOURNAL, MAY 19(14 

be obtained in terms of the perturbed anomaly / if it is understood that 
Oi (l) in the quadrature and in a x is obtained by (28) and if the integrand 
of (31) is multiplied by I /v. 

3.1 Oblateness Effects 

We briefly illustrate some results from Lagrange's method. In the 
well-known example of oblateness perturbations, the first-order solu- 
tions for a, e, and i turn out to be entirely periodic and not very in- 
teresting. 8 The remaining elements, however, exhibit secular terms. 
Using only the ./-term of (9) we find from (21), (22) and (31) 

JR 2 

fix = Qo - —m F^ ( cos *°) Cfi - /o) + periodic terms (32) 

a 2 (l - e o 2 ) 2 



u>! = wo — (cos t ) (fii — fio) + 



JR 2 



a>i 



a 2 (l - eo 2 ) 2 



(l- | sin 8 to) 



(/i - /»> + i n R ,2M X P*- 
a 2 e (l — e 2 ) 2 

n n JR 2 
fix 



(33) 



T\ = To - — ~ if-, 2\3 

ni a (1 — eo j 



X \t{\ + eo cos/) 3 [1-3 sin 2 to sin 2 (a>„ + f)))f\ 
+ — (1 - eo 2 )* [(«i - w) + (cos to)(Oi - Oo)] 



(34) 



JR 2 



2/1 . Z\| 

niao (1 — eo )• 



(l - | sin 2 toj (/, - /o) + p.t. 



where /, /i , /o represent the unperturbed anomaly. ( We have omitted 
the superscript zero for convenience.) Equation (32) confirms the well- 
known secular behavior of the node. It turns out to shift westward for 
< to < tt/2 and eastward for w/2 < to < ir. At t = ir/2 it remains 
stationary, as would be expected from symmetry. 

The secular component of a/ 1 ', according to (33), reduces to the well- 
known term 

JR 2 (5 cos 2 to — 1) 
2a 2 (l - eo 2 ) 2 ' 

It represents an advance of perigee for ^ to < 63°26' and for 116°34' < 
to ^ t- For G3°26' < to < 116°34' perigee regresses, and at the "critical" 
angles 63°26' and 116°34' it is reduced to periodic motions (as far as 
the first-order analysis indicates). We note that the periodic terms in 
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(33) contain e in the denominator, and we expect w to behave un- 
stably for near-circular orbits (as one might expect for geometric 
reasons). Indeed, this singular behavior can be expected also in other 
examples, according to (19), (22) and (31). Furthermore, some diffi- 
culties will arise with small values of i , according to (20) and (21). 
These cases of near-circular and near-equatorial orbits can be accom- 
modated by redefining the orbit elements in various ways. While such 
modified elements are less accessible to a geometric interpretation, they 
do not encumber the calculation of perturbed satellite positions as a 
function of time. For the sake of brevity we must forego additional 
details here. 

3.2 Luni- Solar Gravitation 

We omit a discussion of fix , since it shows periodic perturbations 
only. Substitution of (8) into (19)-(22) and (31) yields 

+ p.t. 

where 

fh = hx p + nhij p + n { z„ 
>k = /j.r 2 + nuy P + n 2 z p 

%i = to + — 3 5 m f° 3 Id + W)hh 2i - (l - owu 

m e r p VI — e 2 

f4(^)-fi};> 

^ = ilo + 3mpao 3 [(l + 4e„ 2 )(Q I - + (1 - e^Qij),] 

2m e r p b sin i (l — eo 2 ) 1 

(37) 



(36) 



-ia^ME + * 



0)1 = 0)o + COS«o(fi — fli) 
Q 3 

+ ^7 (1 - eo)l4/il ~ fh -' p] (38) 



-ife)'-€>- 



862 THE BELL SYSTEM TECHNICAL JOURNAL, MAY 1964 

Tl = To - + (1 ~ ^^ [&i - co„ + (fli - J2o) COS to] 



♦»(-*) 



3Gw p (l - cnT(<i - to) 
2n u r p b (l -f e cos/o) 2 



• [ri - 8(h cos/o + A, Bin/o) 2 ] + ^4 (39) 

AfbQ ill' ]> 

.[a + 3,,= - 3*' (1 + W) - ^ (1 - e„ ! )] 

The subscripts i in (36) and (37) denote partial differentiation with 
respect to the inclination. The secular term of e (1) is a significant feature 
of our present results. It indicates that even near-circular satellite orbits 
can experience an unstable buildup of eccentricity due to luni-solar 
perturbations. The rate of this perturbation is proportional to the factor 
m. p ao 3 /m e r p 3 and is usually very small; moreover, the coordinates x v , 
y p , z p of the perturbing body are really time-dependent, which would 
modify the first-order result. Nevertheless, a long-period change of the 
eccentricity due to luni-solar gravitation has been observed in some 
satellite orbits. 

The explicit form of the periodic terms in a ll) and e (1) contains the 
factor \/e , while O'", to (1) and r (1) contain l/(sin i ). Again this necessi- 
tates the use of specially modified elements for low c and i . One set of 
elements which is particularly suited to the problem of interplanetary 
perturbations is due to Stromgren. He utilized the fact the S2, i, o> are 
nothing but a set of Euler angles orienting a system of orbital coordi- 
nates with one axis through pericenter, one at / = v/2, and one normal 
to the orbit. The rotation of these axes with respect to inertial space 
conveys the same information as the perturbations of S2, i, to. The idea 
is akin to Roberson's method for anticipating secular terms in the 
perturbation of coordinates (Section 4.3 and Ref. 15). 

3.3 Higher-Order Analyses 

The preceding examples are indicative of results to be found in the 
vast literature on perturbations in the osculating elements. We have 
merely covered the gist of this approach and several ideas which will 
be useful in the appraisal of other methods. Some of the better-known 



. (21 



M (2> 
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contributions in terms of osculating elements arc contained in papers 
by Krause, O'Keefe, Kozai, and Iszak. 

In principle the quadratures (18)-(22) and (31) could be evaluated 
iteratively to generate higher-order results. This procedure rests directly 
on Poincare's convergence proof, and a formal technique based on this 
approach is commonly attributed to Poisson. In several aerospace 
publications this has been done to obtain second- and third-order secular 
terms for oblateness effects. The algebraic labor is considerable, though 
typical secular terms such as (40) and (41) tend to be reasonably com- 
pact (see Refs. 7 and 0) : 

( .)ir.J'n Cn /t _ . 2 . • s n \n i • / i n \i- 

^p— ; — : (1-3 sin in sin 0»)[l + e sin (w + 6 { ,)\ 

2a 8 (l - <V) 5 (40) 

X (5 sin 2 i» — 4) cos (w + fti) 

9tJ 2 R 4 sin 2t' (l r i 2 . i/e • 2 • ,-» 

-— -— j^t'o sin 2w„ - |(5 sill In - 4) 

4a 4 (l - c -) 1 

2 • (41) 

[\en sin 2a>(i — <•„ cos a*, (cos 0„ + ^ cos 30o) 

— e u sin a>o (sin U — \ sin 35,,)]). 

These results represent secular increments over a 27r step in 6, the cen- 
tral angle measured from the nolle, where 6 is the initial value of 6. 
Considerable emphasis must be placed in the derivation of such expres- 
sions on checks from the conservation of energy and angular momentum 
and duplicate execution of the algebra. (One likes to think that more 
elaborate explicit expressions will be attainable with the advent of com- 
puter algebra.) A notable contribution in this area was made by Merson, 8 
who presents second-order secular terms for ,/ and first-order secular 
terms for the next four higher harmonics of the earth's potential. He 
also advocates the use of "smoothed" elements which reduce the ampli- 
tude of first-order periodic terms that might otherwise be inimical to 
prediction accuracy. 

IV. PERTURBATIONS IN THE COORDINATES 

We turn now to a description of satellite motions directly in terms 
of the position and velocity vectors. While these are dynamically equiv- 
alent to the instantaneous orbit elements, we note that the time de- 
pendence of the dynamic state variables in this form reflects the anom- 
alistic motion as a primary effect. Therefore the long-time, secular 
changes of the orbit may not be obtainable with the same clarity or pre- 
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cision as in terms of the osculating parameters. On the other hand, the 
position-time history gives a direct account of the satellite motion in 
space and is useful for many aerospace applications. This prospect has 
stimulated several engineering analyses in recent years. 

To the analyst with a general background in mechanics it would seem 
quite natural to approach a system of equations of the type (14) by 
standard perturbation techniques. Thus one could assume the pertur- 
bation series 

,■(/) =x m (t) +« (U (0 +«V*(0 + ••• 

V (() = jj^d) + K// m (/) + «V"(0 + ••• ( 42 ) 

z(l) = z m (t) +kz (, \0 +k 2 /- ,, (0 + ••• 

and determine successively higher-order terms from the appropriate 
governing equations, following essentially Poisson's procedure. The 
traditional Enckc method pursues this line of attack. However, Cartesian 
iuertial coordinates have not enjoyed as much popularity as spherical 
ones, which seem more compliant with the geometry of satellite orbits. 
In the following, therefore, we shall concentrate on reference frames of 
this general type. 

4.1 Perturbations in Equatorial Spherical Coordinates 

A rather well-known perturbation analysis for oblateness effects is 
that due to Blitzer, Weissfeld, and Wheelon. 10 It uses the conventional 
equatorial spherical coordinates, r, <p, \f/ (see Fig. 2) in terms of which 
the equations of motion read: 

r — np' — r cos <p i = — -= — — ~. — li — sin <?J (4.3) 

?- 7'* 

-r O'V) + J ' 2 sin <p cos <pip 2 = ^— sin <p cos <p (44) 

at r 6 

% (r 2 cos%^) =0. (45) 

at 

Here we have considered the first aspherical term in the earth's poten- 
tial only. Since this is a zonal harmonic and does not contain \f/, the last 
equation has a vanishing right-hand side. Then a first integral of this 
equation 

r 2 cos 2 ipyp = p = const., (40) 
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Fig. 2 — Equatorial spherical coordinate system. 

representing the conservation of angular momentum about the polar 
axis, permits a change to \p as independent variable. In addition, it is 
convenient to introduce the definitions 



1/r = u, tan ip = S, and [l/(r cos <p)\ = W. 
Equations (43) and (44) now become 



II'" + w = 



k 3JkR-u 



p 2 (l + S 2 )* 7> 2 (1 + S 2 )* .1 + <S 2 3 



i 



and 



S" + S = 



2.//,/?- Su 



p 2 (1 + , S 2)2 



(47 



(48) 



(49) 



where "primes" denote differentiations with respect to \f/. In this form 
they are readily accessible to a perturbativc procedure. We let 



s = s m + E J n s (n) 

u = >r + f; ./ v 
» =i 

]r m „•"» + g ./"IF"". 



(50) 



Now we note that \p in (46) was understood to represent the actual, 
i.e. perturbed, longitude of the satellite at all times. In order to make 
the connection between this perturbed independent variable and the 
time we find that 
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= t m + E rt M 

i.e. the time itself evolves as a perturbation series. In the zero-order 
solution of (48) to (51 ) the right-hand side of (49) vanishes and we get 

S m = A sin (f + 5), (52) 

where A and 5 are integration constants. On the right-hand side of 
(48) we retain the first term after substitution of S (0) . Then W m follows 
in a straightforward manner. Substituting it into (51), the usual trans- 
cendental expression for time in Keplerian orbits results. It is somewhat 
obscured by the fact that its argument is given in terms of the longitude 
rather than one of the anomalies. We record the simplified form of these 
results for circular orbits, where u = 1/Vo : 

W ™ m I[l+ A 2 sin 2 (* + «)]* (53) 

r 



' - *, tan-'[(l + A 2 ) h tan (* + 8))Y . (54) 

2p ((1 + A 2 )* J*n 

It is important to note that the expressions (52) and (53) represent 
Keplerian (in fact circular) motion only for the unperturbed case: i.e., 
if (54) represents the entire time equation and no higher-order terms 
as in (51) exist. For any perturbed motion, where \p is perturbed in 
relation to time, the zero-order terms of course retain the Keplerian 
forms (52) to (54), but they do not actually represent Keplerian motion. 
If we now proceed to the first-order solution and retain only terms of 
0(./) in (48) and (49), we find 

07 E> 2 C( 0> „, (0) 

C (D" , «(D _ 2/cfl O U /ggN 



and 



W mm + W {x) - " ' 



3/b22V 0)> I 1 - 



+ 



Ll + /! 



(56) 



p 2 [l + 5<°n»Ll + S w 3. 
Since the right-hand side of (55) contains only zero-order quantities, 
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we begin our solution there and the result is 

(I) = MTU sin (<A + 5) , cos (* + 3) 

pV A 2 [l + A- sin 2 (iA + 6)] 1 + A 2 

' l- v/l + 4* tftn ~ 1 [a/1 + ^ tan (,A + 8)] (57) 

sin 2(^ + 8) 



1 + yl 2 sin 2 {$ + 8)J ' 

Here we do not show a complementary solution, since it is of the same 
form as S and can he absorbed with the constants A and 8. We could 
now substitute (57) into (56) to find 11'°' and then use (51) to calcu- 
late the time. However, the inverse tangent in »S " constitutes a secular 
term, which is considered an objectionable feature for some applica- 
tions. 

This disadvantage can be avoided by inserting an additional trans- 
formation between \p and the argument of S { '. Instead of using \p -\- 8 
for the latter let it be 

a = A<A + 5 (58) 

where 

X = 1 + £ ./'X (50) 

and the X„ are constants. This device is commonly attributed to Lind- 
stedt. 12 To obtain the zero-order solution we need only substitute a for 
the angular arguments in (52) and (53). However the equation for »S 
changes significantly, viz.: 



97.7?2r.(0) (0) 

" p a [i + sn 2 



&»' + S w = ,r " ," - 2\ 1 S w " (60) 



where the primes now denote differentiations with respect to a. Thus we 
find 



kR 2 (A 2 cos 2 a - \ - A 2 ) 



S K " = Xi A sin a - — 



pM.r„(l + A 2 )(l + A'-sirta) 

)(61) 

The appearance of the free parameter Xi in this result gives us the 



— coso- 
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opportunity to suppress the secular term. Thus, if we choose 



kR 2 
Xi — 2 /, , A 2\3 (62) 

pr (l + A )- 



(61 ) becomes 



S a) = ^TlTjk { tan_1 K* + A 2 )" tan a] - a}. (63) 

pr (l + A y 

Here we have again omitted all terms of the same form as »S (0> . The net 
contribution from the terms in braces is cyclic and has the period 2ir in a. 
According to (58) and (62) this amounts to a period of 

H>-jwSWI (64) 

in ip. In effect, Lindstedt's transformation distorts the independent 
variable to absorb the secular effect. We shall see more of this later. 

In principle we could transform (56) to a and solve for W in a straight- 
forward manner. However, to simplify the algebra, a redefinition of W 
will be convenient. We may backtrack to the explicit form of (48) in 
terms of S, u, and a. 

[, + S m y u {1) " + 2[1 + S m ']S i0) S m, u™' 

+ [S <0) ' 2 + s (0)2 + lW x) 

3/fefl 2 ,n>2 f 1 2] 

= -tfT u [i + S m * " 3j 

- 2[S (0) 5 (1) + £ (0) '£ (1) ' + 8 m '\]u m 
and take 

JF (1) = (1 + S (0)2 )V o . (66) 

Then we obtain 



(65) 



r(1) „ , W(I) fcflvni - 2^ (o,z ] 



jp^" _|_ n/ tu = 



p 2 A 5 

.W) 



2u [jg (o )jg( i, + ^(.,^(1,/ + g (0)/« XJ 

A 3 



(67) 



where A = (1 + A 2 sin 2 a) . 

Substituting (63) and ignoring the complementary solution for W {1) we 



* Note that the formulas (36) and (38) in Ref . 9 contain several misprints. 
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have 



W ™ = o 3 m TTZXtTi W 4 U' - 1) sin 4 a , RQ . 

3p 3 r 2 (l + A 2 ) 3 A 8 (68) 

+ A 2 (9A 4 + 2/1 2 - 7) sin 2 a + 5A 4 + 2A 2 - 3]. 

Now it only remains to find a relation between a- and time. From (46) 
it is clear that 

t- u> = i (Vcos 2 ^ = i [ —±—-, (69) 

P-Vo P^.w f (l + S 8 )X' v ' 

which we expand up to 0(./). This results in 



'"" " P (i + x')> ten "' 1(1 + Al)> *" '' (70) 



and 



* (,) = 



(71) 



V J, A 2 L A A 2 J 

kR 2 r f(l + A 2 )> . 
- ~ p » (1 + A 2 )* [-1^- ^ sin 2^ + 12(1 + A*)*c] 

+ [2 - 3A 2 + ^ ((1 + .4 2 ) sin 2 a - cos 2 «r)l 

■tan -1 [(1 + A^taiier]} 

and completes this analysis of near-circular orbits. 

Throughout the foregoing discussion we have used a perturbed 
coordinate, i> or <r, as the independent variable. In principle, we could 
have done without Lindstedt's device, and we could have used the un- 
perturbed longitude ^ (0) as the independent variable. This approach 
has the attractive feature that the zero-order solution (in terms of ^ (0) ) 
represents true Keplerian motion. The perturbed longitude could be 
expressed in terms of \f/ (0) as 

* - <r + £, /VW). (72) 

71=1 

However, if one develops the governing differential equations for S (1) 
and W , he discovers that they are completely coupled for this par- 
ticular example. This approach, therefore, loses its practical value. 
In conclusion we note that, in spite of various transformations, the 
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final results (63), (08), (71) of this analysis seem rather awkward, 
considering the fact that they represent the relatively trivial first-order 
oblateness perturbations of a near-circular orbit. This is of course due 
to the choice of spherical equatorial coordinates to represent the motion 
in a nonequatorial orbit. That disadvantage was eliminated in other 
formulations, to be considered next. 

4.2 Perturbations in Orbital Spherical Coordinates 

As the title of this section indicates, it is more natural to use the 
plane of unperturbed motion as the fundamental plane for inclined 
orbits. A typical analysis in this category is that by Anthony, Fosdick 
et al. 13,14 The coordinates r, 0, of their reference frame (see Fig. 3) 
take the place of r, ^, (t/2) - <p in Fig. 2. The angle a = (tt/2) - <p 
(Fig. 3) is introduced occasionally for trigonometric simplifications. 

The left-hand sides of the equations of motion of course are not al- 
tered by this change of coordinates, but the right-hand sides ( represent- 
ing oblateness perturbations) acquire the forms shown in (73) to (75). 
As usual, we introduce u = 1/r and p = r 2 6 = 6/v 2 and change the 
independent variable from t to 0. We note that is the perturbed cen- 
tral angle in the nominal orbit plane. Thus, the equations of motion 
become 

(pit') + pu(/3' 2 + sin 2 0) = (k/p)[\ + ./#V( 1 - 3 cos 2 a)] (73) 

(pd')' - p sin IS cos 13 = {-kJR 2 u/p)[(ain 2 i sin 2 

- cos 2 i) sin 2f3 (74) 

-f- sin 2i cos 2/3 sin 0] 

(p sin 2 )' = ( -A- JR-u/p) [sin 2 i sin 2 sin 20 

(75) 

+ \ sin 2i sin 20 cos 0] 

where primes denote derivatives with respect to 0. We subject, this 
variable to the first-order Lindstedt transformation 

<r = 0(1 + Ai) (76) 

and use the "ansatz" 

u = u (0) (<r) + Ju (i) (a) 

p = p ( V) + W) (77) 

= (tt/2) + J0 O, (<r). 
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Fig. 3 — Orbital spherical coordinate system. 

Without an excessive loss of generality we may assume a horizontal 
launch at the node and, by the definition of the plane as nominal 
orbit plane, the initial velocity vector v lies in it. Thus, at t = 0: r = 
r , f = 0, /3 — ir/2, /3 = 0, and 6 = v u /r . In general, v ( , will be such as 
to produce an elliptic orbit. The zero-order results are 



(0) 

7> = r n /'o 



and 



where 



= (k/r 2 Vo)[i + (' cos a] 



e = (roVo/k) — I 



(78) 



and has the form of a Keplerian eccentricity. As in Section 4.1, the 
zero-order solution (78) will represent Keplerian motion only for the 
unperturbed case, i.e. when a = 6 = 9 . Now the first-order solutions 
follow in a straightforward manner: 

7; (l) = (—k 2 /2r a v 3 )(R/r ) 2 sin 2 i {] -f- fe — c cos a — \e cos 3a 

— COS 2a\ 

u a) = (fc 3 /2 2 /r V){l + he 2 + sin%'(- * + |e - fe 2 ) 

— |(e + sin" *) cos 2a — ^e sin" i cos '3a 

— ■£*<? sin 2 i cos 4o- 

— [I + \e 2 - sin 2 i(l - f<' + |e 2 )]coser}, 



(79) 



80) 
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where we had to select the Lindstedt parameter as 

M = (fc 2 flVr V)(tsin 2 ;- 1) (81) 

to avoid a secular term in (80). Finally 



r2i->2 



(1) = kR sm2i [( ] + } gin a _ a cog a _ i e gin 2ff] (g2) 
2r V 

It is clear that the secular term which was absorbed by the Lindstedt 
parameter has to do with the apsidal precession. In the absence of addi- 
tional Lindstedt parameters we have no countermeasures against the 
secular term in (82), which reflects the nodal regression. (Note that the 
latter was counteracted by the Lindstedt transformation of Section 4.1, 
since it was the only secular effect to be considered for near-circular 

orbits.) 

The time equation for this example can be written in a straightforward 
fashion. From the definition of p it follows that 






" (1 -f x '> a. 



where <t and a x correspond to the time limits / and U ■ An expansion 
to first-order terms yields 



t w + Jl {,) = / 



ffl do- 



(0) (0)2 

o V u 



n i T v a) w (1) l 

- •/ / (0) (0)5 Xi + ~ + 2 —- \da. 



(83) 



This leads to 



foW J — c sin a 



r = 



k 2 \(1 - e 2 )(l + e cosct) 



(84) 



which we recognize as being of strictly Keplerian form but in terms of 
the perturbed angle a. Similarly 



(i) __ R*_ fsin ct(4 — c 2 -f 3c cos a) 
~ Wo \(1 - e 2 ) 2 (l + e cos cr) 2 



«"' = 
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1 + fe + fc 2 + le 3 
+ 1 g- — f + ¥ + ft; - f e I sin « 

2 - h + |e a (85) 



sin a , j| 



(1 - e 2 )(l + e cos a) 



[2(1 -re 

(1 - 



f^'-fey-a:- 



The set of results (78) to (8">) gives a reasonably convenient descrip- 
tion of first-order oblateuess perturbations which might be useful in the 
targeting and guidance of space vehicles. Extensions to near-parabolic 
and hyperbolic trajectories follow quite readily. As in Section 4.1, we 
note that the analysis might have been executed in terms of an unper- 
turbed independent variable, viz. t( " instead of 6, and in that case the 
zero-order solution would represent true Keplerian motion. 

The inclusion of secular perturbations in the independent variable a 
serves the same purpose as the definitions of "mean elements" intro- 
duced by Breakwell et al., by Hansen, in the von Zeipel method, and in 
modern averaging techniques. The Lindstedt transformation is not the 
most powerful device in this category but it can be extended to absorb 
secular effects in more than one coordinate. This will be illustrated in 
the next section in terms of "secular rotations" of the reference frame. 

4.3 Perturbations in Rotating Spherical Orbital Coordinates 

The idea of using suitable coordinate transformations with arbitrary 
parameters to neutralize secular trends was exploited in a more general 
way by R. E. Roberson. 1 His approach uses the orbital coordinates r, 
0, 5 [= (71-/2) — /3], in agreement with Section 4.2, but assumes that 
the entire reference frame will be subjected to three monotonic rota- 
tions, corresponding to three Euler angles, such that the satellite motion 
relative to this reference frame exhibits only periodic perturbations. 
This kinematic outlook on secular trends forms an interesting parallel 
to several classical procedures. Roberson himself makes some illuminat- 
ing comparisons between engineering analyses such as Refs. 9, 10, and 
12 to 1() and traditional formulations in terms of mean variables. He 
restricts his analysis to first-order perturbations, realizing that a con- 
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sistent higher-order theory would have to include contributions from 
other physical effects and various coupling terms. Some of his remarks 
seem quite perspicacious in comparison with the other aerospace litera- 
ture of that time. 

The angular velocities stipulated for the reference frame must of 
course depend on the secular effects that need to be absorbed. In the 
presence of several physical disturbances the different angular motions 
of the coordinate system can be superposed to first order, and the re- 
sultant motion of the reference frame will succeed in neutralizing all the 
secular effects simultaneously. This seems intuitively obvious and can be 
demonstrated in a straightforward fashion. 1 

In Fig. 4 the angles ft and i define a mean orbit plane, in that each of 
them manifests a secular rate. Now the satellite position is given in 
terms of the orthogonal system x, y, 2, which displays a secular variation 
with respect to the node (and this corresponds to the third Eulerian 
rotation). Let the three secular rotations be denoted k$2 (1) , k(iU /dt), 
m w where k is the perturbation parameter. They will in general be 
functions of ft, l and the characteristics of the perturbation source. 

Turning to the problem of first-order oblateness perturbations, we set 
k = J and assume the appropriate form for the perturbing potential. 
Now let 



»(0J 



= 0u + / = 6, + r + 



•// 



id 



(86) 



where / = at t = 0. We adopt / as independent variable and let de- 




Fig. 4 — Rotating spherical orbital coordinates. 
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rivatives with respect to it be denoted by primes. Assuming that the 
secular rates are constants, we let 

q - n (05 + /q ( 7 

l = i m + Ji a) j (87) 

= + t /a>" ) /. 

As usual, the equations of motion arc transformed by means of 

u = 1, /• and p = /•''/ (88) 

and we find: 
(pu'Y + pu[(8' - ST" cos sin if 

+ (cos 8 + u> ( " cos 5 + fi u> cos <p cos pf\ (89) 

- (/,• />)|1 + 3/fiV(l - 3sin>)] = 

\p(8' - Q ll) sin i cos 0)\' + p[( 1 + <S U) ) 2 sin 5 cos 5 

+ (1 + co < ")i'2 <l, (siii 8 cos »> cosy> + cos 5 sin <p) (90) 

+ 12 sin <p cos j/ cos <p] + ( JkR 2 Qu/p) sin ^ cos y> cos j/ = 

[p cos 8 ( 0' cos 5 + D 1 " cos ^ cos v)]' 

+ /;|( I + to'" )<T" sin 5 cos 5 sin i cos 

+ ft ' cos sin / sin 5 cos i» cos y? (91) 

- S'ft (1) sin i'sin 0] 

+ (QJkR'u'p) sin y? cos 5 sin i cos = 0, 

where <p is the latitude anil v is defined in Fig. 5. We have made a slight 
digression from orderly progress in this step by setting i a) = 0. This is 
prompted by previous experience with this problem — viz., that no 
first-order secular perturbations occur in i — and would have developed 
from the later calculations in any event. 
Using the forms 

u = ,r + ju w 

V = //'" + V* (92) 

8 = J8 iu 
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Fig. 5 — Definition of v. 

we reduce (89)-(91) to equations of 0(1) and 0(J). The zero-order 
solution is of course 



s (0) E 0> p (0) = const) 
and 



w ( °> = (fc/p (0)2 )[l +ecos(/- a)]. 



(93) 



As in previous examples, we see that this will represent Keplerian 
motion only in the absence of perturbations, i.e. if / = / and fi = 
£ (1) = 0, yielding an inertial reference frame. We assume that the con- 
stants of integration (p (0) , e, a) are chosen such that (93) with/ = 
yields the satellite position and velocity at t = 0. 

Proceeding with the solutions to 0(J) in the usual fashion, we re- 
quire that 



6 (1) = - 



3Jk 2 R 2 cos i i0) 



V 



(0>« 



(94) 



(95) 



in order to avoid a secular term in 5° and 

a> 0) = (3S 2 fc 2 /2p (0)4 )(5cos 2 t (0) - 1) 

to avoid one in u (1) . These of course reflect the nodal and apsidal pre- 
cessions. The complementary solution for p (1) introduces one constant 
of integration, and the complementary solutions for 5 (1) and u l have 
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the form 

A sin/ +B cos/. (96) 

Since the zero-order solution already accounts for the dynamic state of 
the satellite at t = 0, the first-order solution encounters homogeneous 
initial conditions as far as they do not reflect the rotation of the reference 
frame. Thus at/ = 0: 

u w = 5 ( " = 0, u il) ' = 

5 (n ' - ft (1) sin i m cos do = 0, 

and 

(l// (0) )(2w ( VV 0) + W (0) V ,) ) + a> a) 

+ f2 (,, (cos; (0) - sin i m sin d ) = 0. (9?) 

These govern the first-order constants of integration. [A little reflection 
shows that the forms (96) for 6 (l) and u (1) can be interpreted geometri- 
cally as constant changes of £2 and 6 to 0(J). Roberson anticipates this 
by introducing such constants in (87) and using them in place of two 
of the integration constants for 6 <n and u w . However, nothing seems to 
be gained by this artifice and, if anything, it distracts from a systematic 
procedure.] 

Finally, the time equation follows as usual in terms of/ to 0(J). 
Roberson proceeds to invert it, though the computational gains do not 
seem to justify this algebraic labor. 

So much for our sketch of Roberson's procedure. Its extension to 
higher-order analyses is fairly obvious. At every level of refinement, 
(«/"), three coordinate rotations may be introduced — which are com- 
mensurate with the three degrees of freedom of the satellite problem 
whose secular trends we are trying to neutralize. 

For "medium-range" prediction formulas it seems an open issue 
whether the rationale described in this section and traditional astro- 
nomical devices (like the "auxiliary ellipse" used by Hansen or the von 
Zeipel transformations based on Hamilton-Jacobi techniques) offer a 
computational advantage over the straightforward development of the 
Poisson method for successive higher-order terms. With the advent of 
computer algebra the latter technique may be quite satisfactory for 
many applications. However, for "long-range" predictions and life- 
time studies it seems advisable to employ the accredited astronomical 
techniques of "extracting" secular effects and anticipating long-period 
terms in one way or another. 
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V. MOVING RECTANGULAR ORBITAL COORDINATES 

We close this article with a formulation which calculates the posi- 
tion offsets for a satellite from its unperturbed orbit in an explicit 
form. 18 ' 19 Instead of reckoning the perturbations in terms of the quan- 
tities r, 8, or 5, which are defined relative to the center of the 
earth, we now consider a coordinate system whose origin is the nom- 
inal satellite position 0' on the unperturbed orbit (see Fig. 6). The 
latter may be defined by the initial conditions at time t , viz. r 
and v . We establish an orthogonal triad about the moving point 0' 
with £ in the radial direction, r\ in the direction of anomalistic motion, 
and f normal to the orbit plane. In the guidance engineer's language 




p IG c — Moving rectangular coordinates centered at nominal satellite posi- 
tion. 

these represent offsets in altitude, range, and cross-range. Any non- 
vanishing coordinates in this system are the effects of errors in the 
initial conditions or of geophysical forces. It is clear that this description 
of the perturbed motion can be quite useful in guidance studies, e.g. 
to exhibit the relative motion between a space station (given by 0') 
and a transfer vehicle (located at £, v, f) in a homing maneuver. In the 
subsequent discussion / will always represent the unperturbed true 
anomaly in the nominal orbit and 6 = u + / the unperturbed central 
angle.* No Lindstedt transformations or perturbative coordinate rota- 
tions will be employed to develop this theory into a more sophisticated 
prediction scheme. Instead, we concentrate on the geometric interpre- 
tation of various results. 



* We depart from earlier notations by omitting the superscript (0) from un- 
perturbed quantities for simplicity. 
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The equations of motion can be derived by the standard Lagrangian 
or Hamiltonian formalism, and their linearization to 0(£, tj, f) yields 

f" - 2V - € - 2[f + e(r - V ) sin/J/(l + e cos/) 

q 3 (l - e 2 ) 3 T7 „, 4 (98) 

= r Kf/(] + e cos/) 

V + 2£' - ij - [- v + 2r(7>' + $) sin/]/(l + e cos/) 

= ^- K,/(l + C COS/) 

r + [f - 2f'c sin/].'( 1 + r cos/) 

_ a\l -e 2 )\ 7 .,, , _ .,< (100) 



/,• 



1V(1 + ccos/) 4 



where primes denote derivatives with respect to / and V is the perturba- 
tive potential, which exists in addition to the central body term —k/r. 
The subscripts of V denote partial derivatives with respect to £, y, or f . 
The solution of the homogeneous set (98) and (99), where V m 0, 
represents a complementary solution for the cases where V ^ and 
will be needed to satisfy the initial conditions. For an elliptic nominal 
orbit this solution has the form 



t = 8a 



.1 + 



1 - e 2 3e / k Y , . . ,1 

fHoi/ " 2 1(1 - &W) ° ~ t) Sin/ J 

-5.ac-os/-5r,( 7x A^-)sin/ 

( * v 

\(1 - e*)a) 
o(l - e 2 ) 



(101) 



1 + e cos f y 



(1 + ccos/) + 6o> 



1 + e cos / 



£ = tt^t ^7x [&' si " e ~ ^ sin * cos 0]. (103) 

(1 + e cos/) 

This result can be adapted to hyperbolic, parabolic, and near-parabolic 
orbits without much trouble. The constants of integration 8a ■ • ■ 8t are 
of course just a set of numbers to be determined from the initial condi- 
tions, but the symbols we use for them indicate the parameter changes 
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of the nominal orbit that they represent. Alternatively, these constants 
could be given in terms of £ • ■ ■ to, the perturbations of the position 
and velocity at t . That form is more descriptive for various guidance 
applications. Thus we find for nominally circular orbits 

£ =2W + 4£ - (2W + 3fc) cos (/ - /o) + hi sin (/ - /o) (104) 
V = vo - 2&/ - 3(W + 2fe) (/ - /..) 

+ 2(2W + 3&) sin(/ - /„) + 2fo' cos (/ - / ) 
f = fo ' sin (/ - /,) + fo cos (/ - /o). ( 106) 

(Since a nominal perigee does not exist for this case, we assume that 
the angle w, whose existence is still implied by the notation /, has some 
arbitrary value a. Without loss of generality we can set / = so that 

d = w.) 

In a guidance application So • ■ ■ fo' might represent the errors result- 
ing from a position and velocity determination or a corrective thrust 
maneuver and (104) to (106) would then describe the "run-out" as a 
function of time. In an "orbit sensitivity" study these expressions can 
be used to demonstrate the effect of £ • • ■ fo' on the orbit parameters. 
In a homing maneuver the same expressions would represent the rela- 
tive motion between the two vehicles attempting a rendezvous. In 
principle, two relative position measurements X: £, ij, f at separate 
times suffice to determine all the constants in (104) to (106), and a 
corrective maneuver could be planned to drive the residuals to zero at a 
specified instant or by successive approximations. 

Particular solutions of (98) to (100) can be found in a straightforward 
manner if e = 0. For e^Owe construct these solutions as power series 
in e, for lack of a better expedient. We consider the series to 0(e) and let 
them be denoted by 

i 

l-li + * + e£u ./= 1,2,3 (107) 

3 

f = fl + f2T«S I'J 
3 

where 

\fi£x = the complementary solution (104) to (106) 
£>i?2?2 = a particular solution representing V to 0(k) 
fiQihi = solution reflecting e-(h, Vi , h) on the right-hand sides of 
(98) to (100) 
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faith = solution reflecting e-(& , m , f 2 ) on the right-hand sides 
/:if/:i/';t = solution reflecting <'•(!'{, V v , l' f ) on the right-hand sides. 
The following explicit general forms can be given for these solutions: 



V* 



h = (a/lc) 1-2 J f, df + 2 cos/ / \\ cos f df 

+ COS / J ?i sin / r// + 2 sin / J P, sin / df 

- sin/ / I'{ cos / df 

= (a 3 /k) h |j P, df df + 2 / f { d/ - 4 sin/ | f , cos/ df ( 108) 

- 2 sin / J l\ sin / df + 4 cos / J V v sin / df 

- 2 cos/ / P { cos fdf 



1 1 >«. 1 



f, ■ („/•//,) Los/ / P r sin / # - sin / f P r cos / # 

where we take /o = for the lower limit of all quadratures. 
/, = ( v „ - 2£„') sin / - $(V + 2| ) cos/ - 3(W + 2&) 

•'•sin/ 
f/i = 7(W + 2fr) sin/ + (ij„ - 2t,/) cos/ - 3(W + 2$ 

•>' cos/ - (6,'/2) cos 2/ - ( V + #&>) sin 2/ 

fei = - (fo/2) - (fn'/2) sin 2/ - (ft 2) cos 2/. 

The terms /-. , (/■• , and h-> are obtainable from expressions analogous to 
(108) but with the following substitutions: 

2(£/ - m) sil1 / - -h cos/ for ( -a/k)~\\ 

2(^' + fc ) sin / -|- ^j cos / for ( -a 1 /*) P, (110) 

2f/ sin / + f, cos / for ( - a/k) f f 

and /shafts follow from (108) if we substitute 

(4a 8 A ) f s ,„. f for ( -a 3 /fc) f { ,„, f . (Ill) 

Since the differential operators for all of these partial solutions are of 
the form 
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f - 2 V ' - 3€ 

i* + 2? (ii2) 

r + r, 

no explicit complementary solution of 0(e) is provided: i.e., the con- 
stants £» • • • £</ will be used to satisfy the i.c.'s to all levels of accuracy. 
These will differ from zero only if the nominal orbit is taken to differ 
from r and v at i . 

Specific results may now be obtained by the above formulas, which 
lend themselves to a geometric interpretation of perturbed satellite 
motions. For the oblateness effect one finds 

fe = (JR 2 /a)[-l + sin 2 i(f + * cos 20)] 

ij 2 = (JR 2 /a)[(2 - 3 sin 2 i)f + tV sin 2 i sin 26] (113) 

f s = (JR 2 /2a) sin 2i[/cos B - \ sin 0]. 

The terms in 26 reflect the doubly symmetric distortion of the orbit 
due to the oblateness of the gravitational field. The constant term in 
|a and the secular term in ij-2 reflect the additional mass of the equatorial 
bulge. Combining (113) with (104) to (IOC)) into a complete solution, 
we observe that the constant term in £ is 

Art - (JR 2 /a)[-l + | sin 2 i\ + 2i? ' + 4&, 

and the secular term in -q (114) 

o,A0 = f[(JR 2 /a)(2 - 3 sin 2 i) - 3(W + 2&)], 

which represent the differences between the nominal circular orbit 
and the mean circular orbit resulting from the perturbations. Since 
Art = and Ad = do not yield linearly independent conditions for 
£ and t? ', we cannot effect a launch so that the radius and the mean 
angular rate coincide with the nominal ones (determined for a spherical 
earth) unless sin i = y/2~JZ. On the other hand, it turns out that we 
can preserve the nominal inclination of the orbit by choosing f = and 

fo' = (JR 2 /2a) sin 2% cos O . (115) 

Now, if we designate A? 2 = [fdSlo', we find for the nodal regression 

ft = - A '^ k \ = - n j f^Y C os i (116) 

sin 1 2x0 s \ a / 

and this agrees with the well-known result. 
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In the case of drag perturbations one replaces T~ ? , V v , Vj of (110) 
by the appropriate components of (11): 

F $ = 

F n = - (CV4p /2m)[(/,-/a)" - <m cos i\(k/a) h (117) 

F; = — (C D Ap /2m) a sin i(ka)' cos 6 

and finds 

f 2 = (2a 3 /k)FJ 

ij, = (a 3 F„// C )[4 - f/ 2 ] (118) 

f, = (a 3 F r /4/.-)[2/ tan - cos 0]. 

Noting that 

1 r^/,9=2 I r A" C nAfwrfff Sill » /„m 

~ — If 2 J 9=0 — -,t = — -. — n ' Uiy; 

2ira df 4mk' 

we have agreement with standard results for the orbital precession 
due to diurnal winds. 

If we extend this drag analysis to 0(e) we find 

/» + /:. = (aF v /2k)\% sin/ - ./' cos/ - :if sin/] 

o 2 + o 3 = (o'F-/2A;)[6/ sin /" + 1) cos f - Sf cos f] 

' _ (120) 

fh + /':i = [aFf/{4k cos 0)][:> cos (to + 2/ ) — V co - s w — ,/ sin d> 

-/sin (to + 2/)], 

which are simple enough to permit a further extension to cases where 
po = p(/) is variable around the orbit. The details are straight forward.' 

It is of course understood that any of these results should be ac- 
companied by lirjifi if a general solution is desired. This, however, 
adds nothing to the characteristics of a particular perturbation. The 
formulas (104) to (100) and (108) to (111) can also be applied to a 
variety of other effects such as luni-solar perturbations and radiation 
pressure. 

The motivation behind the results of this section was to give a geo- 
metrically tangible account of perturbed satellite motions over a frac- 
tional period or just a few periods. This may be useful in various tar- 
geting, intercept, and rendezvous operations. On the other hand, the 
formulations of Sections III, 4.2, and 4.3 form the beginnings of ephem- 
eris computing techniques and orbit lifetime studies. These subjects 
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have been pursued further in several higher-level methods (see Section 
I), some of which deal partly with the elements and partly with co- 
ordinates and make occasional use of contact transformations. They 
may be considered a stepping stone to full-fledged astronomical pertur- 
bation analyses, about which there exists an extensive literature. 
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